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i-^ ■ Abstract. We use the Hough transform to analyze data from the second science 

run of the LIGO interferometers, to look for gravitational waves from isolated 
pulsars. We search over the whole sky and over a large range of frequencies 
and spin-down parameters. Our search method is based on the Hough transform, 
which is a semi-coherent, computationally efficient, and robust pattern recognition 
technique. We also present a validation of the search pipeline using hardware 
signal injections. 
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1. Introduction 



This paper presents some partial results for a wide parameter space search for periodic 
gravitational waves using data from the LIGO detectors. The most promising sources 
for such waves are isolated pulsars. Previous searches for gravitational waves from 



^Jy pulsars have been of two kinds. The first is a search targeting pulsars whose parameters 

are known through radio observations. These searches typically use matched filtering 
techniques and are not very computationally expensive. An example of such a search 
is PP which targets PSR J1939+2134 using data from the first science runs of the 
' LIGO and GEO detectors. The end result is an upper limit on the strength of 

the gravitational wave emitted by this pulsar and therefore on its ellipticity. See 
also |2] which applies some of the techniques presented in pQ to a large number of 
known pulsars using data from the second science run of the LIGO detectors. The 
second kind of searches look for pulsars which have not yet been observed by radio 
telescopes. This involves searching over large parameter space volumes and turns 
out to be computationally limited. This is because looking for weak continuous wave 
signals requires large observation times to build up signal to noise ratio and to claim a 
detection with some degree of confidence; on the other hand, the number of templates 
that must be searched over, and therefore the computational requirements, increase 
rapidly with the observation time. An example of such a search is |3j where a 2-day 
long data stretch from the Explorer bar detector is used to perform an all sky search 
in a narrow frequency band around the resonant frequency of the detector. 

All the searches mentioned above rely on a coherent integration over the full 
observation time; it is well known that this is the optimal method. However, a full 
coherent integration is computationally expensive and it is therefore also useful to 
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consider methods which are less sensitive but computationally inexpensive. Such 
methods typically involve semi-coherent combinations of the signal power in short 
stretches of data. The Hough transform is an example of such a method Using 
this method, we perform an all-sky search over a large frequency range using two 
months of data from the LIGO detectors. As in all the searches mentioned above, we 
assume that the pulsar does not glitch during the full observation time considered. 

Section [21 briefly describes the waveforms that we are looking for. Section 
describes our search method, the Hough transform. The search pipeline and the 
parameter space we search over is given in section 0] The search results are given 
in section [5j Section presents a validation of our search method using hardware 
injected signals and finally section [7| concludes with a summary of our results and 
plans for further work. 



2. The expected waveform 

The form of the gravitational wave emitted by an isolated pulsar, as seen by a 
gravitational wave detector is [Sj : 

h(t) = F + (t,i>)h + (t) + F x (t,ip)h x (t) (1) 

where t is time in the detector frame, ip is the polarization angle of the wave, F + _ x 
are the detector antenna pattern functions for the two polarizations. If we assume 
the emission mechanism is due to deviations of the pulsar's shape from perfect axial 
symmetry, then the gravitational waves are emitted at a frequency which is twice the 
rotational rate f r of the pulsar. Under this assumption, the waveforms for the two 
polarizations h+ tX are given by: 

1 ~\~ COS 2 L 

h + = ho cos <!>(£), h x = ho cos tsin $(t) , (2) 

where i is the angle between the pulsar's spin axis and the direction of propagation of 
the waves, and ho is the amplitude: 

16tt 2 G I zz ef? 

ho~— —■ (3) 

Here d is the distance of the star from Earth, I zz is the z-z component of the star's 
moment of inertia with the z-axis being its spin axis, and e is the equatorial ellipticity 
of the star. The phase <&(i) takes its simplest form in the Solar System Barycenter 
(SSB) frame where it can be expanded in a Taylor series. Up to second order: 

$(t) = $ + 2tt (j (T - T ) + |/(T - Tof^j . (4) 

Here T is time in the SSB frame and To is a fiducial start time. The frequency fo 
and the spin-down parameter / are defined at this fiducial start time. In this paper, 
we include only one spin-down parameter in our search, i.e. we ignore the higher 
order terms in Eq. This is reasonable because, as we shall see below in section^] 
our frequency resolution is too coarse for the higher spindowns to have any effect for 
reasonable values of the pulsar spindown age. 

Neglecting relativistic effects which do not affect us significantly in this case (again 
because of the coarseness of our frequency resolution), the instantaneous frequency 
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f(t) of the wave as observed by the detector is given, to a very good approximation, 
by the familiar non-relativistic Doppler formula: 

/W-/»=/» X ^ (5) 

c 

where t is time in the detector frame, v(t) is the velocity of the detector at time t, n 
is the direction to the pulsar, f(t) is the instantaneous signal frequency at time t and 
is given by 

fit) = fo + f it - 1 ) . (6) 

Equations JSJ) and 10 describe the time-frequency pattern produced by a signal, and 
this is the pattern that the Hough transform is used to look for. 



3. The Hough transform 

The Hough transform was invented by Paul Hough in 1959 as a method for finding 
patterns in bubble chamber pictures from CERN [Hj and it was later patented by IBM 
0. The Hough transform is also well known in the literature on pattern recognition 
to be a robust method for detecting straight lines, circles etc. in digital images; see 
e.g. [Sj for a review in this field. A detailed discussion of the Hough transform as 
applied to the search for continuous gravitational waves can be found in 0j . A closely 
related semi-coherent method is the stack-slide algorithm described in [Sj. 

The idea of the Hough transform can be illustrated by the following simple 
example. Consider the problem of trying to detect straight lines in a noisy two- 
dimensional digital image. The digital image is assumed to be made up of pixels 
which can be in one of only two possible states, namely "on" or "off". Let (x, y) be 
the coordinates of the center of a typical pixel. We are looking for a pattern which is 
parameterized by two numbers (m, c) such that 

y = rax + c. (7) 

The parameter space (ro, c) is assumed to be suitably digitized so that it is also made 
up of pixels. To find the most likely value of (m, c), we proceed as follows. For each 
pixel (£, y) which is "on" , we mark all the possible values of (m, c) which are consistent 
with it, i.e. we mark all pixels in the (m, c) plane lying on the straight line y = mx + c 
with a This is repeated for every pixel which is "on". The end result is an 

integer, the number count, for every pixel in the (m, c) plane. In the case when the 
digital image is too noisy and no straight lines can be detected, the number counts 
would be uniformly distributed in the (m, c) place. The presence of a sufficiently 
strong signal would lead to a large number count in at least one of the (m, c) pixels, 
and the largest number count would indicate the most likely parameter space values. 

This method enables us to mark all the possible templates consistent with a 
given observation without stepping through the parameter space point-by-point. This 
leads to a significant gain in computational speed. Furthermore, each observation, no 
matter how noisy, only adds at the most +1 to the final number count. These two 
features are the chief virtues of the Hough transform method: computational speed 
and robustness. On the other hand, the Hough search is likely to be less sensitive 
than the stack-slide search considered in [Jj]. The tradeoffs between sensitivity versus 
efficiency and robustness are yet to be studied in detail and will be important in the 
context of a hierarchical search |3 ^] . 
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In our case, the Hough transform is used to find a signal whose frequency evolution 
fits the pattern produced by the Doppler shift © and the spin-down ||BJ| in the time- 
frequency plane. The parameters which determine this pattern are (/ ,/,n); a point 
in this four dimensional parameter space will be denoted by £. This parameter space 
is covered by a discrete cubic grid whose resolution is described in section The 
result of the Hough transform is a histogram, i.e. an integer (the number count) 
for each point of this grid. The starting point for the Hough transform are N short 
stretches of Fourier transformed data; each short stretch will be called an SFT (Short 
Fourier Transform). Each of these SFTs is "digitized" by setting a threshold p th on 
the normalized power pi~ in the k th frequency bin:. 

Pk — TP, a It \ ■ \°> 

1 cah^n(Jk) 

Here Xk is the value of the Fourier transform in the k th frequency bin corresponding 
to a frequency fk, T coh is the time baseline of the SFT, and S n (fk) is the single-sided 
power spectral density of the detector noise at the frequency fk- We require that T coh 
is small enough so that the signal does not shift by more than, say, half a frequency 
bin within this time duration. For frequencies of ~ 300Hz, this restricts T coh to be 
lesser than ~ 60min 0]. In this paper, we work with SFTs for which T coh = 1800s. In 
principle, we could choose T coh to be greater, but we are restricted by the duty cycle 
of the interferometers in that we should be able to find suitably long time periods 
during which the detector is in lock. Furthermore, the data should be stationary over 
the chosen time period. The choice of 1800s is a suitable compromise for all the three 
interferometers during the S2 run. 

This thresholding produces a set of zeros and ones (called a "pcakgram" ) from 
each SFT. This set of pcakgrams is the analog of the digitized two-dimensional image 
described earlier. The Hough transform is used to calculate the number count n at 
each parameter space point starting from this collection of peakgrams. Let p(n) be 
the probability distribution of n in the absence of a signal, and p(n\h) the distribution 
in the presence of a signal h(t). It is clear that < n < N, where N is the number 
of SFTs, and it can be shown that for stationary Gaussian noise, p(n) is a binomial 
distribution with mean Nq where q = e~ Pth is the probability that any frequency bin 
is selected: 

P(n)= ( N n ^jq n (l-q) N - n . (9) 

In the presence of a signal, the distribution is ideally also a binomial but with a slightly 
larger mean TV rj where, for weak signals, r\ is given by 

/ ? = 9 {l + ^A + 0(A 2 )} . (10) 

A is the signal to noise ratio within a single SFT, and for the case when there is no 
mismatch between the signal and the template: 

with h(f) being the Fourier transform of the signal h(t) (see 0] for details). The 
approximation that the distribution in the presence of a signal is binomial breaks 
down for reasonably strong signals. This happens mainly due to two reasons: i) the 
random mismatch between the signal and the template used to calculate the number 
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count and ii) the amplitude modulation of the signal which causes r\ to vary from one 
SFT to another and for different sky locations. The result of these two effects is to 
"smear" out the binomial distribution in the presence of a signal. 

Candidates in parameter space are selected by setting a threshold n th on the 
number count. The false alarm and false dismissal rates for this threshold are defined 
respectively in the usual way: 

N n th -l 

a= P= £ p H^)- (12) 



n=0 



We choose the thresholds (n th ,p th ) based on the Neyman-Pearson criterion of 
minimizing (3 for a given value of a. It can be shown 0] that this criterion leads, 
in the case of weak signals (i.e. A << 1), large N, and Gaussian stationary noise, 
to p th rs 1.6. This corresponds q = e~ Pth ks 0.20, i.e. we select about 20% of the 
frequency bins from each SFT; for weak signals, this turns out to be independent of 
the choice of a and signal strength. Furthermore, n th is also independent of the signal 
strength and is given by: 

nth = Nq + ^2Nq(l-q) crfc" 1 (2a) (13) 

where, as before, q = e~ Pth and erfc -1 is the inverse of the complementary error 
function. These values of the thresholds lead to a certain value of the false dismissal 
rate f3 which is given in 0] . The value of /3 of course depends on the signal strength, 
and it turns out that the weakest signal which will cross the above thresholds at a 
false alarm rate of 1% and a false dismissal rate of 10% is given by 



_ 8.54 S n (f ) , . 

Equation l|14|) gives the smallest signal which can be detected by the search, and is 
therefore a measure of the sensitivity of the search. 

The data analyzed in this paper correspond to LIGO's second Science Run (S2) 
that was held for 59 days, from February 14 to April 14, 2003. The GEO detector 
was not running at that time, but all three LIGO detectors were operating with a 
significantly better sensitivity than during the first science run. The LIGO detectors 
comprise one 4 km facility in Livingston, Louisiana, (LI) and two, 4 km and 2 km 
respectively in Hanford, Washington (HI and H2); see eg. [TT]. For our purposes, we 
note that the duty cycles of the detectors during the S2 run were 37% for LI, 74% for 
HI and 58% for H2. The number N of 30min SFTs available for LI data were 687, 
1761 for HI and 1384 for H2. 

Figure ^ shows the expected sensitivity for the Hough search by the three LIGO 
interferometers during the S2 run. Those ho values correspond to the amplitudes 
detectable from a generic source with a 1% false alarm rate and 10% false dismissal 
rate, as given by Eq. I|14|l. It should be kept in mind that Eq. I|14J) significantly over- 
estimates the sensitivity of the search for unknown pulsars because it does not include 
the mismatch between the signal and the template. Furthermore, due to the large 
number of templates involved in the search, a false alarm rate of 1% is too large in 
practice, and it would result in too many potential candidates. A false alarm rate of 
~ 10~ 13 would be more realistic since this would lead, in the ideal case, to less than 
one candidate over the parameter space points considered in this search. 

Assuming that the gravitational wave emission mechanism is due to deviations of 
the pulsar's shape from perfect axial symmetry, from Eq. ©. (Eq. (2.9) in 0j) and 
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Figure 1. Typical sensitivities of the three LIGO detectors during the S2 run 
with a 1% false alarm rate and 10% false dismissal rate. 



Eq. 1)14(1 . we can estimate the nominal astro-physical reach of the search for the three 
detectors: 



d 



m^GN^I^efi T coh 



8.54c 4 
10 45 g-cm 2 , e 



(15) 



For a value of I zz = 10 45 g-cm 2 , e = 1CT 5 , and for typical parameters of the S2 run, 
this corresponds to a distance of about 20-30 parsecs. It should be kept in mind 
that this is not a realistic figure for the astrophysical reach of the search; it does not 
consider the mismatch between the template and signal, and it does not use the more 
realistic false alarm rate mentioned above. Due to these effects, it turns out that 
Eq. ((ToT) overestimates the astrophysical reach by a factor of about 2-3. 



4. The search pipeline 

Data from each of the three LIGO interferometers are used to analyze the same 
parameter space region. This section describes the portion of the parameter space 
(/o,/,n) that we search over, and the resolution of our grid in this portion of the 
parameter space. 

The total observation time is approximately T ob3 ps 5.2 x 10 6 sec corresponding to 
the S2 science run. We search for pulsar signals in the frequency range 200-400Hz with 
a frequency resolution: 5f = T coh ~ = 5.556 x 10 _4 Hz. The resolution Sf in the space 
of first spindown parameters is given by the smallest value of / for which the intrinsic 
signal frequency does not drift by more than a single frequency bin during the total 
observation time: Sf = 5f ■ T~j ~ 1.1 x 10 _10 Hz/s. We choose the range of values 
—/max < / < 0, where / max = 1.1 x 10 _9 Hz/s. This yields eleven spin-down values. All 
known pulsars (except for a few supernova remnants) have spindown parameters less 
than this value. This value of / max is equivalent to looking for pulsars whose spindown 
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Figure 2. Pipeline for the semi-coherent Hough search for a single interferometer. 



age t = f I f is at least 1.15 x 10 4 yr. This also shows that the approximation to drop 
higher spindown terms in Eq. (0} is reasonable; with a spindown age of 1.15 x 10 4 yr 
as above, we would need a total observation time of ~ 10 yr for the second spindown 
to cause a frequency drift of half a frequency bin. 

The resolution 89 in sky positions is frequency dependent, with the number of 
templates increasing with frequency and is given by 89 = ^(<5#) min , where (<5#) min is 
given m Eq. (4.14) of g|. This yields a resolution of about 9.3 x 10 _3 rad at 300Hz. 
This resolution corresponds to ~ 1.5 x 10 5 sky locations for the whole sky. 

The pipeline used to search over this parameter space is described in figure [3 
The figure is divided into four distinct blocks. The top-left block is the preparation 
of the SFTs: the data stream is broken up into segments, calibrated, and a discrete 
Fourier transform is applied to each segment. The calibration connects the error-signal 
from the interferometer to the actual value of the strain, and this is calculated in the 
frequency domain. These SFTs are passed onto an optional conditioning step. This 
is meant to remove any known spectral disturbances from the SFTs. In the present 
paper, we only present results for which no data conditioning is applied to the SFTs. 
Finally, data from the three interferometers are analysed separately. 

The rest of the pipeline consists of two conceptually distinct parts: the actual 
Hough search and the process of setting upper limits. The Hough search has been 
described earlier; a threshold is set on the normalized power of each SFT, replacing 
thereby the SFTs by a set of pcakgrams. In this paper, we only present partial results 
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Figure 3. Top: maximum, mean, minimum and standard deviation of the 
number count of all the Hough maps in the frequency band 206-207 Hz. The 
data corresponds to LI for the entire S2 run using 687 SFTs with a time baseline 
of 30 minutes. Bottom: the solid line corresponds to the LI number-count 
distribution obtained in that band, and in red circles the theoretical expected 
binomial distribution for 687 SFTs and a peak selection probability of 20%. 



from this Hough search and not the process of setting upper limits in any detail, except 
to say that this is the conventional frcqucntist upper limit based using Monte-Carlo 
simulations. We set upper limits in each 1Hz frequency band, based on the loudest 
event observed in that band. This will be presented elsewhere. 

5. Partial results from the search code 

As described earlier, the first step in this semi-coherent Hough search is to select 
frequency bins from the SFTs by setting a threshold on the normalized power defined 
in Eq. JHJ. This requires a reliable estimate of the power spectral density S n for 
each SFT, for which we employ a running median applied to the periodogram of each 
individual SFT. The running median is a robust method to estimate the noise floor 
|12| which has the virtue of discarding outliers which appear in a small number of bins, 
thereby providing an accurate estimate of the noise floor in the presence of spectral 
disturbances and possible signals. 

As an illustrative example, some results of the Hough search in a 1 Hz frequency 
band arc shown in figure^ The first panel of this figure shows, for every frequency bin, 
the maximum, minimum, mean and standard deviation of the number counts for all 
sky-locations and all spindown values. As expected, the mean is approximately Nq = 
0.2 x 687 « 137. Similarly, as expected, the standard deviation is y q(l — q)N « 10. 
The second panel of figure shows the distribution of number counts in this band and 
compares it with the expected binomial distribution in the absence of any signal. We 
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Figure 4. Graph of the LI maximum number count per frequency analyzed, 
maximized over all spin-down values and sky locations. The dash-dotted line is 
the corresponding threshold n th for a false alarm a of 10~ 10 . 



find excellent agreement with the expected binomial distribution, and this is true in 
all frequency bands which are relatively free of spectral disturbances. 

Figure 01 shows the largest number count obtained for every frequency bin (i.e. 
the maximum number count over all sky-locations and spindown values for a given 
frequency value). As this figure shows, several environmental and instrumental noise 
sources are present. The sources of these disturbances are mostly understood. They 
consist of broad 60Hz power lines, multiples of 16Hz due to the data acquisition system, 
and the violin modes of the mirror suspensions in a neighborhood of 345Hz. The 60Hz 
lines are rather broad, with a width of about ±0.5Hz, while the 16Hz data acquisition 
lines are confined to a single frequency bin. In addition to the above disturbances, we 
also observe a large number of multiples of 0.25Hz. While these lines are known to be 
instrumental, their exact physical origin is yet to be determined. 

6. Pipeline validation with hardware signal injections 

Two artificial pulsar signals were injected for a duration of 12 hours at the end of the 
S2 run into all three LIGO interferometers. These injections were designed to give 
an end-to-end validation of the search pipeline starting from as far up the observing 
chain as possible. 

The two artificial signals were injected at frequencies of 1279.123Hz (PI) and 
1288.901Hz (P2) with spindown rates of zero and — 10 _8 Hzs _1 respectively, and 
amplitudes ho of 2 x 10 -21 . The signals were modulated and Dopplcr shifted to 
simulate sources at fixed positions on the sky with ip = 0, cost = and <fi = 0. PI 
was injected at a right acsension of 5.147 rad and a declination of 0.3767 rad, while 
P2 had a right ascension of 2.3457 rad and a declination of 1.2346 rad. 

The resolution in the space of sky positions and frequencies are the same as in 
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Figure 5. Hough maps for the hardware injected signal PI in LI. Map 2442 (top- 
left) corresponds to 1279.123333 Hz, and contains the template which is closest 
to the signal. The top-right panel is a zoom of this map, showing the signal 
more clearly. Maps 2222 and 2662 (bottom-left and bottom-right) have a larger 
mismatch in frequency; they correspond to 1279.112222 Hz and 1279.134444 Hz 
respectively. The signal is detected in these maps also, but with a mismatched 
sky-location. PI was injected at a right acsension and declination of 5.147 rad 
0.3767 rad respectively. This sky-location corresponds roughly to the center of 
the skypatches shown in these figures. 



section 01 but the spin-down resolution depends on the total observation time, and 
this now turns out to be -2.28624 x 10~ 8 Hzs" 1 for LI, -1.77024 x 10~ 8 Hzs" 1 for 
HI, and —1.93533 x 10~ 8 Hzs -1 for H2. As before, for each intrinsic frequency we 
analyze 10 different spin-down values. The portion of the sky analyzed has a width of 
0.5 radians x 0.5 radians centered around the location of the injected signals. 

Figure shows some results for pulsar PI using LI data. There were 14 SFTs 
available in the duration when the pulsar was injected. The top left and top right 
panels of figure shows Hough maps corresponding to the frequency and spin-down 
values nearest to the injected signal. Although the presence of the signal is clearly 
visible, 12 hours of observation time is not enough to identify the location of the source 
in the sky. In particular, while the signal is identified with a high number count in 
these hough maps, one can still identify the signal in Hough maps corresponding to 
different frequencies and spindowns with high number-counts, but with a mismatch 
in the sky location. This is shown in the bottom left and bottom right panels of [S] 
Similar results were found for pulsar P2 and the other detectors, thus providing an 
important validation of this search pipeline. 
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7. Conclusions 

In this paper, we have described the idea of the Hough search and the search pipeline 
used to analyze data from the second science run of the LIGO interferometers. We 
have shown some outputs of the Hough search pipeline in the frequency range 200- 
400Hz, over the whole sky, and the first spindown parameter. We have also validated 
the search pipeline by showing that the search can detect hardware injected pulsar 
signals. Work is in progress to compute astrophysical upper limits using the search 
pipeline presented in this paper. 

The eventual role of the Hough transform is in a hierarchical scheme HOj . 
The Hough transform could be used as a computationally inexpensive and robust 
method for quickly scanning large parameter space volumes and producing significant 
candidates for a follow-up search using a more sensitive method. 
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